home *** CD-ROM | disk | FTP | other *** search
/ Super PC 34 / Super PC 34 (Shareware).iso / spc / UTIL / DJGPP2 / V2 / DJTST200.ZIP / tests / libc / ansi / math / elefunt / machar.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-08-23  |  10.8 KB  |  428 lines
  1. #include "elefunt.h"
  2.  
  3. #ifdef __STDC__
  4. void store(float* p);
  5. #else
  6. void store();
  7. #endif
  8.  
  9. void
  10. machar(ibeta, it, irnd, ngrd, machep, negep, iexp, minexp, maxexp,
  11.     eps, epsneg, xmin, xmax)
  12. int *ibeta,
  13.     *iexp,
  14.     *irnd,
  15.     *it,
  16.     *machep,
  17.     *maxexp,
  18.     *minexp,
  19.     *negep,
  20.     *ngrd;
  21. float *eps,
  22.     *epsneg,
  23.     *xmax,
  24.     *xmin;
  25. /***********************************************************************
  26. #
  27. #     THIS SUBROUTINE IS INTENDED TO DETERMINE THE CHARACTERISTICS
  28. #     OF THE FLOATING-POINT ARITHMETIC SYSTEM THAT ARE SPECIFIED
  29. #     BELOW.  THE FIRST THREE ARE DETERMINED ACCORDING TO AN
  30. #     ALGORITHM DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951,
  31. #     INCORPORATING SOME, BUT NOT ALL, OF THE IMPROVEMENTS
  32. #     SUGGESTED BY M. GENTLEMAN AND S. MAROVICH, CACM 17 (1974),
  33. #     PP. 276-277.  THE VERSION GIVEN HERE IS FOR SINGLE PRECISION.
  34. #     CARDS CONTAINING  CD  IN COLUMNS 1 AND 2 CAN BE USED TO
  35. #     CONVERT THE SUBROUTINE TO DOUBLE PRECISION BY REPLACING
  36. #     EXISTING CARDS IN THE OBVIOUS MANNER.
  37. #
  38. #
  39. #       IBETA   - THE RADIX OF THE FLOATING-POINT REPRESENTATION
  40. #       IT      - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING-POINT
  41. #                 SIGNIFICAND
  42. #       IRND    - 0 IF FLOATING-POINT ADDITION CHOPS,
  43. #                 1 IF FLOATING-POINT ADDITION ROUNDS
  44. #       NGRD    - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION.  IT IS
  45. #                 0 IF  IRND=1, OR IF  IRND=0  AND ONLY  IT  BASE  IBETA
  46. #                   DIGITS PARTICIPATE IN THE POST NORMALIZATION SHIFT
  47. #                   OF THE FLOATING-POINT SIGNIFICAND IN MULTIPLICATION
  48. #                 1 IF  IRND=0  AND MORE THAN  IT  BASE  IBETA  DIGITS
  49. #                   PARTICIPATE IN THE POST NORMALIZATION SHIFT OF THE
  50. #                   FLOATING-POINT SIGNIFICAND IN MULTIPLICATION
  51. #       MACHEP  - THE LARGEST NEGATIVE INTEGER SUCH THAT
  52. #                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, EXCEPT THAT
  53. #                 MACHEP IS BOUNDED BELOW BY  -(IT+3)
  54. #       NEGEPS  - THE LARGEST NEGATIVE INTEGER SUCH THAT
  55. #                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, EXCEPT THAT
  56. #                 NEGEPS IS BOUNDED BELOW BY  -(IT+3)
  57. #       IEXP    - THE NUMBER OF BITS (DECIMAL PLACES IF IBETA = 10)
  58. #                 RESERVED FOR THE REPRESENTATION OF THE EXPONENT
  59. #                 (INCLUDING THE BIAS OR SIGN) OF A FLOATING-POINT
  60. #                 NUMBER
  61. #       MINEXP  - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT
  62. #                 FLOAT(IBETA)**MINEXP IS A POSITIVE FLOATING-POINT
  63. #                 NUMBER
  64. #       MAXEXP  - THE LARGEST POSITIVE INTEGER EXPONENT FOR A FINITE
  65. #                 FLOATING-POINT NUMBER
  66. #       EPS     - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH
  67. #                 THAT  1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER
  68. #                 IBETA = 2  OR  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
  69. #                 OTHERWISE,  EPS = (FLOAT(IBETA)**MACHEP)/2
  70. #       EPSNEG  - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT
  71. #                 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2
  72. #                 OR  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
  73. #                 OTHERWISE,  EPSNEG = (IBETA**NEGEPS)/2.  BECAUSE
  74. #                 NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT
  75. #                 BE THE SMALLEST NUMBER WHICH CAN ALTER 1.0 BY
  76. #                 SUBTRACTION.
  77. #       XMIN    - THE SMALLEST NON-VANISHING FLOATING-POINT POWER OF THE
  78. #                 RADIX.  IN PARTICULAR,  XMIN = FLOAT(IBETA)**MINEXP
  79. #       XMAX    - THE LARGEST FINITE FLOATING-POINT NUMBER.  IN
  80. #                 PARTICULAR   XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
  81. #                 NOTE - ON SOME MACHINES  XMAX  WILL BE ONLY THE
  82. #                 SECOND, OR PERHAPS THIRD, LARGEST NUMBER, BEING
  83. #                 TOO SMALL BY 1 OR 2 UNITS IN THE LAST DIGIT OF
  84. #                 THE SIGNIFICAND.
  85. #
  86. #     LATEST REVISION - OCTOBER 22, 1979
  87. #
  88. #     AUTHOR - W. J. CODY
  89. #              ARGONNE NATIONAL LABORATORY
  90. #
  91. ***********************************************************************/
  92. {
  93.     int i,
  94.         iz,
  95.         j,
  96.         k,
  97.         mx;
  98.     float a,
  99.         b,
  100.         beta,
  101.         betain,
  102.         betam1,
  103.         y,
  104.         z;
  105.     float t1,
  106.         t2,
  107.         t3;            /* temporaries for subexpressions */
  108.  
  109. #ifdef MSC
  110.     /* Need to allow for overflow here--MSC otherwise aborts job */
  111.     _control87(CW_DEFAULT | EM_OVERFLOW, 0xffff);
  112. #endif
  113.  
  114. #ifdef __TURBOC__
  115.    /* Get around DOS 3.2 bug that crashes PC after 8 floating-point
  116.    traps, as described in README file for Turbo C 1.5 distribution.
  117.    This call disables all floating-point traps. */
  118.    _control87(MCW_EM, MCW_EM);
  119. #if 0
  120.    _control87(IC_AFFINE,IC_AFFINE);    /* Turbo C 1.5 uses IC_PROJECTIVE */
  121.                     /* Microsoft C 5.0 uses IC_AFFINE */
  122. #endif
  123. #endif
  124.  
  125.  
  126.     /* determine *ibeta,beta ala Malcolm */
  127.  
  128.     a = ONE;
  129.     do
  130.     {
  131.     a = a + a;
  132.     t1 = a + ONE;
  133.     t2 = t1 - a;
  134.     t3 = t2 - ONE;
  135. #ifdef DEBUG_MACHAR
  136.     printf("L1: %15.7e %15.7e %15.7e %15.7e\n",a,t1,t2,t3);
  137. #endif
  138.     }
  139.     while (t3 == ZERO);
  140.  
  141.     b = ONE;
  142.     do
  143.     {
  144.     b = b + b;
  145.         t1 = a + b;
  146.     t2 = t1 - a;
  147. #ifdef DEBUG_MACHAR
  148.     printf("L2: %15.7e %15.7e %15.7e\n",b,t1,t2);
  149. #endif
  150.     }
  151.     while (t2 == ZERO);
  152.  
  153.     t1 = a + b;
  154.     *ibeta = INT(t1 - a);
  155.     beta = FLOAT(*ibeta);
  156.  
  157.     /* determine *it, *irnd */
  158.     *it = 0;
  159.     b = ONE;
  160.     do
  161.     {
  162.     (*it)++;
  163.     b = b * beta;
  164.     t1 = b + ONE;
  165.     t2 = t1 - b;
  166.     t3 = t2 - ONE;
  167. #ifdef DEBUG_MACHAR
  168.     printf("L3: %15.7e %15.7e %15.7e %15.7e\n",b,t1,t2,t3);
  169. #endif
  170.     }
  171.     while (t3 == ZERO);
  172.  
  173.     *irnd = 0;
  174.     betam1 = beta - ONE;
  175.     t1 = a + betam1;
  176.     t2 = t1 - a;
  177.     if (t2 != ZERO)
  178.     *irnd = 1;
  179.  
  180.     /* determine *negep, *epsneg */
  181.  
  182.     *negep = *it + 3;
  183.     betain = ONE / beta;
  184.     a = ONE;
  185.  
  186.     for (i = 1; i <= *negep; ++i)
  187.     a = a * betain;
  188.  
  189.     b = a;
  190.     for (;;)
  191.     {
  192.     t1 = ONE - a;
  193.         store(&t1);
  194.     t2 = t1 - ONE;
  195.         store(&t2);
  196.     if (t2 != ZERO)
  197.         break;
  198.     a = a * beta;
  199.         store(&a);
  200.     (*negep)--;
  201. #ifdef DEBUG_MACHAR
  202.     printf("L4: %15.7e %15.7e %15.7e %15.7e\n",a,t1,t2,t3);
  203. #endif
  204.     }
  205.  
  206.     *negep = -*negep;
  207.     *epsneg = a;
  208.     if (!((*ibeta == 2) || (*irnd == 0)))
  209.     {
  210.     a = (a * (ONE + a)) / (ONE + ONE);
  211.         store(&a);
  212.     t1 = ONE - a;
  213.         store(&t1);
  214.     t2 = t1 - ONE;
  215.         store(&t2);
  216.     if (t2 != ZERO)
  217.         *epsneg = a;
  218.     }
  219.  
  220.     /* determine *machep, *eps */
  221.  
  222.     *machep = -*it - 3;
  223.     a = b;
  224.     for (;;)
  225.     {
  226.     t1 = ONE + a;
  227.         store(&t1);
  228.     t2 = t1 - ONE;
  229.         store(&t2);
  230.     if (t2 != ZERO)
  231.         break;
  232.     a = a * beta;
  233.         store(&a);
  234.     (*machep)++;
  235. #ifdef DEBUG_MACHAR
  236.     printf("L5: %15.7e %15.7e %15.7e %15.7e\n",a,t1,t2);
  237. #endif
  238.     }
  239.     *eps = a;
  240.     if (!((*ibeta == 2) || (*irnd == 0)))
  241.     {
  242.     a = (a * (ONE + a)) / (ONE + ONE);
  243.         store(&a);
  244.     t1 = ONE + a;
  245.         store(&t1);
  246.     t2 = t1 - ONE;
  247.         store(&t2);
  248.     if (t2 != ZERO)
  249.         *eps = a;
  250.     }
  251.  
  252.     /* determine *ngrd */
  253.     *ngrd = 0;
  254.     t1 = ONE + *eps;
  255.     store(&t1);
  256.     t2 = t1 * ONE;
  257.     store(&t2);
  258.     t3 = t2 - ONE;
  259.     store(&t3);
  260.     if ((*irnd == 0) && (t3 != ZERO))
  261.     *ngrd = 1;
  262.  
  263.     /* Determine *iexp, *minexp, *xmin.
  264.        Loop to determine largest i and k = 2**i such that (1/beta) ** (2**(i))
  265.        does not underflow. Exit from loop is signaled by an underflow. */
  266.  
  267.     i = 0;
  268.     k = 1;
  269.     z = betain;
  270.     for (;;)
  271.     {
  272.     y = z;
  273.     store(&y);
  274.     z = y * y;
  275.         store(&z);
  276.  
  277.     /* check for underflow here */
  278.     a = z * ONE;
  279.         store(&a);
  280.     t1 = a + a;
  281.     store(&t1);
  282.     if ((t1 == ZERO) || (ABS(z) >= y))
  283.         break;
  284.     i++;
  285.     k = k + k;
  286. #ifdef DEBUG_MACHAR
  287.     printf("L6: %15.7e %15.7e %15.7e %15.7e %d\n",a,t1,y,z,k);
  288. #endif
  289.     }
  290.     if (*ibeta != 10)
  291.     {
  292.     *iexp = i + 1;
  293.     mx = k + k;
  294.     }
  295.     else
  296.     {                /* for decimal machines only */
  297.     *iexp = 2;
  298.     iz = *ibeta;
  299.     while (k >= iz)
  300.     {
  301.         iz = iz * *ibeta;
  302.         (*iexp)++;
  303.     }
  304.     mx = iz + iz - 1;
  305.     }
  306.     for (;;)
  307.     {    /* loop to determine *minexp, *xmin
  308.            Exit from loop is signaled by an underflow. */
  309.     *xmin = y;
  310.     y = y * betain;
  311.         store(&y);
  312.  
  313.     /* check for underflow here */
  314.     a = y * ONE;
  315.         store(&a);
  316.     if (((a + a) == ZERO) || (ABS(y) >= *xmin))
  317.         break;
  318.     k++;
  319.  
  320. #ifdef DEBUG_MACHAR
  321.     printf("L7: %15.7e %15.7e %d\n",a,y,k);
  322. #endif
  323.  
  324. #ifdef __TURBOC__
  325.     /* Turbo C 1.5 computes minexp = -1074, but the run-time
  326.     library turns 2**p (p < -1023) into 0.0, so we reset the
  327.     limit so as to prevent the use of denormalized numbers.
  328.     Microsoft C 5.0 correctly generates denormalized numbers down to
  329.     2**(-1074). */
  330.     if (k >= 1023)
  331.         break;
  332. #endif
  333.     }
  334.     *minexp = -k;
  335.  
  336.     /* determine *maxexp, *xmax */
  337.     if (!((mx > (k + k - 3)) || (*ibeta == 10)))
  338.     {
  339.     mx = mx + mx;
  340.     (*iexp)++;
  341.     }
  342.     *maxexp = mx + *minexp;
  343.     /* adjust for machines with implicit leading bit in binary significand */
  344.     /* and machines with radix point at extreme right of significand */
  345.  
  346.     i = *maxexp + *minexp;
  347.     if ((*ibeta == 2) && (i == 0))
  348.     *maxexp = *maxexp - 1;
  349.     if (i > 20)
  350.     *maxexp = *maxexp - 1;
  351.     if (a != y)
  352.     *maxexp = *maxexp - 2;
  353.  
  354. #ifdef DEBUG_MACHAR
  355.     printf("L7a: minexp = %d\tmaxexp = %d\n",*minexp,*maxexp);
  356. #endif
  357.  
  358.     *xmax = ONE - *epsneg;
  359.     store(xmax);
  360.     if (*xmax * ONE != *xmax)
  361.     *xmax = ONE - beta * *epsneg;
  362.     store(xmax);
  363.     *xmax = *xmax / (beta * beta * beta * *xmin);
  364.     store(xmax);
  365.     i = *maxexp + *minexp + 3;
  366.     for (j = 1; j <= i; ++j)
  367.     {
  368.         if (*ibeta == 2)
  369.         *xmax = *xmax + *xmax;
  370.         else
  371.         *xmax = *xmax * beta;
  372.     store(xmax);
  373. #ifdef DEBUG_MACHAR
  374. #ifdef IEEE
  375.         /* See remark below */
  376. #else
  377.     printf("L8: %15.7e %d\n",*xmax,j);
  378. #endif
  379. #endif
  380.     }
  381.  
  382. #ifdef IEEE
  383.     /* On a host with IEEE arithmetic, the above code computes 
  384.        *minexp = -1074 (sp: -149) (with gradual underflow; else
  385.        -1023 (sp: -126)), which is acceptable, but it finds
  386.        *maxexp = 3021 (sp: 362?) (with gradual underflow),
  387.        instead of 1024 (sp: 128?)  because the code fragment
  388.  
  389.        if (!((mx > (k + k - 3)) || (*ibeta == 10)))
  390.        {
  391.       mx = mx + mx;
  392.       (*iexp)++;
  393.        }
  394.  
  395.        doubles mx from 2048 to 4096 (sp: 256 to 512).  The
  396.        following code simply compares z with (z*beta)/beta, with
  397.        z = beta**k, until it finds a mismatch, at which point we
  398.        have found the real *maxexp. */
  399.        
  400.     z = beta;
  401.     *maxexp = 1;
  402.     for (;;)
  403.     {
  404.     if (z != z)
  405.     {
  406.       printf("?ERROR: z*beta has generated a NaN instead of infinity!\n");
  407.       break;
  408.     }
  409.     y = z * beta;
  410.     if (y == z)
  411.         break;
  412.     if ((y / beta) != z)
  413.         break;
  414.     (*maxexp)++;
  415.     z = y;
  416. #ifdef DEBUG_MACHAR
  417.     printf("L9: %15.7e %d\n",z,*maxexp);
  418. #endif
  419.     }
  420.     z *= (ONE - *epsneg);
  421.     store(&z);
  422.     z *= beta;
  423.     store(&z);
  424.     *xmax = z;
  425. #endif
  426.     return;
  427. }
  428.